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Abstract 

We present a method to derive the analytical expression of the roughness of a fractal surface 
whose dynamics is ruled by cellular automata. Starting from the automata, we write down the the 
time derivative of the height's average and variance. By assuming the equiprobability of the surface 
configurations and taking the limit of large substrates we find the roughness as a function of time. 
As expected, the function behaves as ir when t <C i x and saturate at w s when t^>t x . We apply 
the methodology to describe the etching model Q|, however, the value of j3 we obtained are not 
the one predicted by the KPZ equation and observed in numerical experiments. That divergence 
may be due to the equiprobability assumption. We redefine the roughness with an exponent 
that compensate the nonuniform probability generated by the celular automata, resulting in an 
expression that perfectly matches the experimental results. 
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I. INTRODUCTION 



The study of stochastic process has accelerated in the last decades, connected to disci- 
plines such as economy, biology, meteorology and neuroscience. Much of the recent advances 
are due to the availability of fast and affordable computer clusters. Such development has 
became feasible, for instance, numerical simulation of large particle systems obeying simple 
repetitive rules which mimetize complex systems. Emergent information and properties have 
been obtained through several techniques and approaches. 

The surface growth phenomena, when treated as a stochastic process, encompass a wide 



application field. Some examples of growth systems are corrosion 



fire propagation 



lu 



ar automata 



6|, lZ|, analytical 



4|, atomic deposition [5|, evolution of bacterial colony nd ce 

models js|. Models have been proposed and studied through experiments [l, 
calculations [9|, and computational simulations [l, 2\. 

In this work we obtain the roughness evolution analytically, for systems of dimension 
1+1, referring to directions perpendicular and parallel to the substrate. Far from being 



mathematical idealizations, these systems have physical meaning. 



'henomena with that 



dimensionality include bacterial colony growth on Petri dish |lONl2|. paper burning, ink 
diffusion on paper and turbulence of liquid crystals 3, 1^]. In the case of paper, por 
example, the burned or stained frontier of the paper can be represented by hf(x,t) where 
< x < L and L is the sheet width, being x and measured along the directions mentioned 
above. The subscript / specify the reference frame fixed in a corner of the sheet. 

Surfaces with different internal dynamics lead to distinct profiles, which can be charac- 
terized by different measures, the most important being the mean value and the standard 
deviation of the surface height. When related to surfaces, the standard deviation is often 
called roughness, defined as 



w(L,t) 




[hf(x,t) -h(t)] 2 dx. 



(1) 



Even if h(t) increases continuously due to the growth process, the dynamic equilibrium 
lead to roughness saturation after a period of roughening buildup. The saturated roughness 
often is a function of the substrate size as the power law w s ~ L a , a being the saturation 
exponent. The saturation occurs at a characteristic time (t x ), and follows the power law 
t x ~ L z , where z is the dynamic exponent. Before saturation (t <C t x ), w(L,t) evolves as a 



power law with the growth exponent 0, w(L, t) ~ t& These properties were incorporated 
in the Family- Vicsek scaling relation 15] 

w(L,t)cxL a f(^), (2a) 



where 



x 13 when i«l 

(2b) 

const when x ^> 1 

Scaling techniques applied to the growth equations may be used to find the aforemen- 
tioned exponents for certain universality classes and dimensions. Growth systems with the 
same exponents are considered to belong to the same universality class, connecting systems 
that seem unrelated. Some well known universality classes are the Edwards- Wilkinson and 



the KPZ 



9]. 



II. METHOD FOR OBTAINING THE ROUGHNESS EQUATION 

In this section we present a method to obtain the equation describing the time evolution 
of the roughness. In the following section we will apply it to the etching model, notwith- 
standing, this method is quite general and can be extended to other models. 



A. The evolution of the mean squared roughness 

The first step is to compute the change in the squared roughness, defined as w q (t) = w 2 (t), 
when the deposition occur in the site % of a substrate with roughness w. The squared 
roughness is equal to mean value of hi 2 (t), with hi defined at the reference frame of the 
mean height 

h t (t) = h/(t)-h(t), (3) 

where % is the particle position. 

Our methodology assume that the change will depend only on the nearest neighbors of site 
i, therefore the variation of the squared roughness will be written as Aw q (w, /it-i, hi, hi + i). 
That assumption defines the celular automatas it can be applied to. We will do the average 
over the ensemble of all possible configurations with roughness w, defining p(w, faj-i, hi, h i+ i) 



3 



as the probability of the values hi-i, hi and hi+i for a given value of w. The evolution of 
the mean squared roughness is 



\ At I J-yJl^fi J-y/Lw 2 -h i+1 2 J-y/Lw2-h l+1 2 -hi 2 At 

p(w, hi_i, hi, h i+l ) dhi_i dhi dh i+1 . (4) 

The integration limits encompass the configurations allowed by the definition of roughness, 
eq. ([1]). As will be discussed latter, the celular automata may prevent the occurrence of 
some configurations, which must have their probabilities assigned to zero. 



B. The rate of change of the quadratic roughness 



Our approach involves calculating the increment of roughness when one iteration is per- 
formed. We use a discrete substrate with sites of length u = 1 and evaluated the system 
before and after the deposition of one particle of height Ay = 1. The squared roughness is 
affected by the following changes of hi 

-i 2 



with 



Before: hi 2 (t) = h{(t) - h(t) 

r 1 2 

After: h?(t + At) = h{ (t + At) - h(t + At) 

= [hi f (t) + Ahi(t)-h(t) -Ah(t)\ 
= hi 2 (t) + Qi. 

Qi =2hi(t)Ahi(t) -2hi(t)Ah(t) 

- 2Ahi(t)Ah(t) + [A^(t)] 2 + [Ah{t)] 2 . 

Using the definition of roughness squared, w q (t) = w 2 (t), we write 



Before: w q (t) = -jr^/i; 2 (t). 

i 

After: w q {t + At) = - ^ [h 2 {t) + Ql ] 



(5) 



(6) 



(7) 
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The unity of time we use corresponds to L iterations, therefore, At = 1 / L and the rate of 
change of w q {t) is 

= -L[Ah(t)} 2 + I £ {^mKit) + [Ah^t)} 2 } , (8) 

i 

where we have used ^ /ij(t) = and Ahi(t) = LAh(t). 

Eq. OH]) is a general formula for the increment of quadratic roughness, independent of the 
iterative algorithm. In order to obtain the roughness for a specific algorithm, it is necessary 
to know the values of Ahi(t) and Ah(t). 

On these grounds, each model results in different values for Aw q (w, hi-!, h iy h i+ i) / At and 
p(w, hi-i, hi, hi+i), which must be deducted for each case. In the next two subsections we 
will assume the equiprobability of the accessible configurations to eliminate the dependence 
on p(w, hi-t, hi, h i+ i). 



C. The equiprobability of the configurations 

Before proposing an expression p(w, hi-!, hi, hi + i) we must remember that there is a finite 
number of possible substrate configuration for each value of w. This finite number is the 
result of the restrictions of hi imposed by eqs. (ED) and (J3J): 

h 2 + h 2 2 + • • • + h L 2 =Lw 2 (9a) 
h + h 2 + ■ • ■ + h L =0 (9b) 

Eqs. ([9]) define a hyperplane and the surface of a hypersphere of radius wyfL both in L- 
D, i.e., in a space of L dimensions. From the intersection of these two subspaces, a spheric 
surface results, which is (L — 2)-D. For each combination of hi-!, hi, and faj+i, the remaining 
L — 3 h's form another spheric surface, now with dimension L — 5. The area of these spheric 
surfaces, with L — 2 and L — 5 dimensions, will be called, respectively, At and A p . 

While all possible surface configurations of a given value of w belong to the (L — 2)-D 
surface, the subset of them for which the values of the triad hi-!, hi, and h i+1 are known 
belongs to the (L — 5)-D surface. We will assume the equiprobability of the configurations 
allowed by eqs. (jHJ), consequently, the probability of a given triad is 

A 

p(w, hi-!, hi, h i+1 ) = (10) 
5 



The assumption of the equiprobability of the configurations, used when deriving that equa- 
tion, disregard two important properties of the dynamics, which is not harmless. The first 
is the inaccessibility, by the celular automata, of several configurations allowed by eqs. (JSJ). 
The second is the different probability of each configuration generated by the celular au- 
tomata. 

Eq. fflOl) can be rewritten by substituting the expression for the area of a hypersphere, 



p(w, hi-i, hi, h 



i+l) 



C D L— 5 
JL-5 tip 

OL-2 tiT 



(11) 



where Sjy are constants with depend on the dimension D, and R p and Rt are the radius of 
the corresponding hyperspheres. To make the notation less clumsy, we will assume, in the 
remaining of this subsection, a deposition at i — 2, implying that only the sites i — 1, 2 ou 
3 may be affected. 

The plane of eq. (I9b|) contains the center of the sphere of eq. fl9a|) . therefore, the sphere 
defined by their intersection also has radius Rt = w\/L. The value of R p may be obtained 
by rewriting eqs. ([9]) as 



K 2 + h 2 + ... + h L 2 =Lw 2 -h\-h\-h\ 



h d + h* 



hi 



(hi + h 2 + h 3 



(12a) 
(12b) 



i.e., a hypersphere with superficial area A p of dimension L — 5. These two equation may be 
combined as 

L 

) 2 , (13) 



- 2 ^ h i h j = Lw<2 - h \- h l ~ h\- {hi + h 2 + h s 

i,j=4 

which we rewrite in a matricial form: 



h/i h 5 



-1 
-1 



1 



hi 
h 5 

hr. 



The eigenvalues of the above square matrix are A 4 = — (L — 4), A 5 = 1, A 6 = 1 , . . . , 
Xl = 1. From these eigenvalues we can define a linear transformation to the set of variables 



h' A , . . . , h' L that eliminated the crossed terms, 



(L - A)h' 2 + h' 2 



h' 2 = Lw 2 - hi 2 - h 2 2 - h 3 2 - (hi + h 2 + h 3 ) 2 . (15) 



Transformed varible h' 4 is equal to the left hand side of eq. (I12bl) 

1 



h' 



1 



With this transformation, eq. (TH1) becomes 

,/ 2 r 2 i. 2 



(/i4 + h + ... + h L ) 

(h 1 + h 2 + h 3 ). 



,/2 ,/2 r 2 t2 t2 t 2 (hl + h 2 + h 3 ) 2 
h 5 + ■ h/l L = LKJ — All — fi 2 — ^3 



L -3 



The radius of this hypersphere is given by 



Rr. 



Lw 2 - h t 2 - h 2 2 - h 3 



{h + h 2 + h 3 ) 
L-3 



21 



(16) 



(17) 



(18) 



We can now rewrite eq. ( TTTj) with the values R? and R p in the asymptotic case L — > oo 

[Lw 2 - h 2 - h 2 2 - h 3 2 ] ^ 



p(w,h 1 ,h 2 ,h 3 ) = rj(L) 
with t)(L) = S L _ 5 /S L _ 2 . 



(Lw 2 



(19) 



D. The roughness evolution with the equiprobability assumption 

A more convenient expression for roughness squared, eq. (J4]), is possible by doing the 
change of coordinates 

hi-i = \/~Lw sin p cos 9 

hi = \Jl~w sin p sin 9 cos tp ■ (20) 
hi+i = \J~Lw sin p sin 9 sin <p 
with the variables defined in the intervals 



7T 



0<p<-, 0<9<n, < <p < 2vr. 
The probability (Tl9|) expressed in these coordinates is 

p(iu, p, 9, tp) = rj(L) (Lw 2 )" 1 cos L ~ 5 p. 



(21) 



(22) 



We can now rewrite the evolution of the squared roughness averaged over the ensembles, eq. 
(H: 



r 2 * r ft A Wq (w, P ,e,<p) . 2 L _ 4 . 

r/(L) / / / — t sin p cos p sin 6 1 dp rfb 1 a<p (23) 

At / Jo Jo Jo At 



'a^\_^ m r 27r r r* Aw q (w, P ,e,<p) n . 

'o Jo Jo 

In the above expression we employed the jacobian 



d/ij d/ij+i = (-^w^ 2 ) 2 sin 2 p cos p sin 6 1 dp d9 dip. (24) 

While eq. @ is exact, eq. (123]) is an approximation based in the equiprobability assump- 
tion. Since we eliminated the dependence on p(w, hi-i, hi, h i+ i), the left hand side of that 
simplified equation depends on the roughness only through the term 

Aw q (w,p,6,<p) 

At ' (25) 

facilitating the resolution of the roughness equation. The celular automata used in each 
model will determine the value of this term, calculate through eq. (jHJ). 



III. EXAMPLE: ETCHING MODEL 

We will know calculate the evolution of the mean value of w q , eq. ( l23l . for the etching 
model. The etching model proposed in 2001 by Mello, Chaves, and de Oliveira [lj] belongs 
to the KPZ universality class 9| . The model mimics the corrosion of a crystal surface by a 
solvent. A time evolution At is defined as one iteration of the following celular automata: 



1. Randomly choose a site i G [1..L]. 

2. If hUit) < h{(t) do hU(t + At) = h{(t). 

3. If h{ +l (t) < h{(t) do hlS + At) = h{(t). 

4. Do h{(t + At) = h{(t) + l. 



The algorithm implements a cell removal probability that is proportional to the number 
of the exposed faces of the cell, a reasonable approximation of the etching process. It could 
also describe a deposition where each exposed face has the same attachment probability. 
For that reason, it can be referred either as particle removal or deposition. 

The scaling exponents found by Mello et al. in 1 + 1 dimension were a = 0.491 and 
(3 = 0.330, placing the model within the KPZ universality class (a = 1/2, (3 = 1/3) [l|. 
Other studies analyzed the model in 1 + 1 and 2 + 1 dimensions, focusing on aspects such 
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FIG. 1: Effects of deposition of a cell in site 2, classified in four cases. 



as dynamic behavior of the roughness and comparisons with other models belonging to the 
KPZ universality class ^,Q-2Q|. 

Before applying eq. (jHJ) to the etching model, it is necessary to know the values of 
Ahi(t) and Ah(t) for one iteration. While the first may be directly obtained from the model 
rules, the second must be separated in four possibilities, depending on the neighbors of the 
deposition site i. These four situations are shown in Fig. (CD) together with the added cells 
when site % is selected. The effect of them on Ah(t) are 



Ah(t) 



Case 1: -^Ay. 

Case 2: \ [Ay + % - h^)] . 



Case 3: \ [Ay 
Case 4: } [Ay 



(hi - h i+1 )] . 

(hi - hi-!) + (hi - h i+1 )] 



(26) 



The rate of change of w q , eq. flS}, can be separated in four parts corresponding to the 
cases of eq. ff26l) . We represent each of them as ^^^-j , with j = [1..4]. From eq. (120]) we 

conclude that hi ~ w, so, eq. (jHJ) implies that ^^ £ j is a quadratic function of w. Since 



the integral (T23]) is a function of w only through Aw q /At, all dependence on w is within 
that term, resulting in a quadratic equation for the roughness dynamics, 

' AW„ 



At 



-c a w 



c b w - c c , 



(27) 



where the minus signals were include to simplify forthcoming calculations. In the limit 
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At = 1/L — > we can do 



At / dt dt 
By replacing the Eq. (|27j) in this expression, we obtain 



A Wq \^dw 1 = 2w dw 



2w dw 2w dw A dw B dw , , 

-c a dt= — — = -. = + , (29) 

w z + — w + — (w — Wi)(w — W 2 ) W — Wi W — W 2 

where W\ and w 2 are the roots of w 2 + + ^ = 0, and A = 2wi , B = ~ 2w ' 2 . The 

solution of the above diferencial equation is 



Wq — WiJ \Wq-W 2/ 

being wq the initial roughness. 

We have found an implicit equation for the roughness as a function of time. Knowing 
that that equation and the Family- Vicsek scaling have both the characteristic time t x , we 
conclude that t x = l/c a . We will show briefly that wiw 2 < 0, and, if we choose W\ > 0, 
then w 2 < 0, A > and B > 0. The limit w — > w s when t — > 00 establishes Wi = w s . 
Finally, roughness w depends on L through the constants t x and w s , therefore, if we want 
the remaining of the equation to be independent of L, we must have w 2 = —Xw s being A a 
positive proportionality factor independent of L. Incorporating that reasoning in eq. ( |30l ). 



and considering a flat initial substrate, we can write 

2A 

e"* /4x . (31) 



w \ L + x ( 1 w \ 



w s ) \ Xw s 

The initial growth, i.e., when w <C w s , is controlled by the exponent 0. In order to make 
explicit that dependence, we expand the left hand side of Eq. (131 p up to the second order, 
which result in 

l + i^4 + 0(3) = e-'/'x. (32) 

The inversion of that equation is 

w(t) &w s y/\(l -e^) 1/2 . (33) 

This expression for w is real only if A > 0, implying Wiw 2 < 0, as we said. From the limit 
of the above expression, 

/ t \ 1/2 

hmw(t) = w a y/X — , (34) 
t— >o \t x J 

we find 13 = 1/2. 

It is interesting to make it clear that if A = 1, Eq. (133]) is identical to Eq. ( 13~T|) . not an 
approximation of it. 
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FIG. 2: Two possible configurations that can not be generated by the etching model. 



A. The aftermath of equiprobability assumption 

The exponent — 1/2 found in the previous section is not the one expected for the KPZ 
university class, which we know the etching model belongs to. The main approximation 
to be blamed for that disagreement is the equiprobability assumption, explicitly, that all 
configurations not forbidden by Eqs. (Q are allowed, and that they have all the same 
probability. 

In growth models based in cellular automata, the internal rules of the model forbid certain 
configurations. For instance, the algorithm that governs the etching model can never lead 
to the configurations shown in the Fig. According to the the model, when one particle 
is deposited on the top of the second site, the first and third sites grow up to the earlier 
height of the second site, which is not the case of Fig. [2j It results that these are prohibited 
configurations. 

Beside that, there is no reason to suppose that the dynamics results in equal probabil- 
ity for the allowed configurations. Indeed, numerical experiments with the etching model 
demonstrated that p(w, hi, hi+i) doesn't agree with eq. (fl9|) . 

The equiprobability assumption disregard the ban of some configurations and the 
non-uniformity of the probability distribution, allowing us to express the probability 
p(w, hi-i, hi, h i+ i) as the ratio of the area of the partial hypersphere (defined by w, hi-i, hi, 
hi+i) and the area of the total hypersphere (defined by w). However, the resulting expression 
for w, Eq. [33j nas n °t the expected exponent (3. 

Aiming to circumvent the effect of the equiprobability approximations, we write the real 
roughness of the etching model, w e , as a function of the the roughness w, which does not 
include those features, 




(35) 
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where the superscript refers to saturation. There are good reasons to reduce all the possi- 
bilities to this form. First, w must increase monotonously with w e , they must agree for null 
roughness, and they have their maximum at saturation. This is the simplest form that that 
fulfil those requirements. Second, for a hypersphere of Euclidean dimension d, and radius 
~ R, the surface area grows as R d ~ l , while surfaces such as those originated by etching, or 
domains as in phase transition are expected to grow as Sj oc , where df is the fractal 

dimension. I.e., they must have a fractal character and dimension df. In this way, Eq. ( )35|) 
is not only the simplest form, it is the only possibility. 
If we substitute eq. (13"o"j) in eq. ( 13T?j) we can write 

where /3 — ^. If we explicitly make w e s oc L a and t x oc L z , that expression becomes one 
possible form of Family- Vicseck scaling relation, Eq. fl2]). 

Although that equation is an expansion strictly valid only for t -C t x it nicely fit to 
the results of the numerical simulations in the whole range of t, as can be seen in Fig. 
03]). Each curves of that figure was obtained by fitting its parameters (3, t x , and w s to the 
corresponding data points. Data from simulation of other surface evolution models (RSOS, 
Edwards- Wilkinson) and higher dimensions have shown the same exceptional agreement 
with Eq. (J36J. 

IV. CONCLUSIONS 

In this work, we presented a method to obtain the roughness equation of the models. 
The method is based on the ratio of the hypersphere areas, interpreted as the probabil- 
ity of occurrence of determined configurations of the interface. The hyperspheres method 
has potential of application in automata cellular models which depend only on the nearest 
neighbors but it needs to be built differently for each type of algorithm. The algorithms 
need to act only in the nearest neighbors and the systems need to be one-dimensional. If 
approximations similar to equiprobability assumption are necessary to solve these models, 
transformation similar to Eq. ( 13 5 p may be necessary. 

The equation possesses three parameters, each of them is connected with one of the 
growth exponents. The modified equation is relevant not only because it fits well to the 



1 — exp ( — 



t 



(36) 
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FIG. 3: Roughness as a function of the time. The points are the results of numerical simulation of 
etching model for several substrat lenghts. The curves are fittings of Eq. (|36|) to the data points. 

data, but also because we can also be used as fitting function for obtaining the value of 
the main parameters of the model, w e s , /3, and t x . In correlated stochastic phenomena , 
analytical results are rather difficult to obtain. In this way we hope that this work may 
inspire research into those systems where even not exact solutions can be considered major 



results (21 1 . 
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